JB-Scent-Net

knitr::opts_chunk$set(echo = TRUE, comment="  ", cache=TRUE, warning=FALSE, message=FALSE, fig.width=14, fig.asp=.678)
library(openxlsx)
library(data.table)
library(RColorBrewer)
pal <- brewer.pal(5,"Set2")
#library(heatmaply)
library(bipartite)

## Loading required package: vegan

## Loading required package: permute

## Loading required package: lattice

## This is vegan 2.6-4

## Loading required package: sna

## Loading required package: statnet.common

## 
## Attaching package: 'statnet.common'

## The following objects are masked from 'package:base':
## 
##     attr, order

## Loading required package: network

## 
## 'network' 1.18.2 (2023-12-04), part of the Statnet Project
## * 'news(package="network")' for changes since last version
## * 'citation("network")' for citation information
## * 'https://statnet.org' for help, support, and other information

## sna: Tools for Social Network Analysis
## Version 2.7-2 created on 2023-12-05.
## copyright (c) 2005, Carter T. Butts, University of California-Irvine
##  For citation information, type citation("sna").
##  Type help(package="sna") to get started.

##  This is bipartite 2.19.
##  For latest changes see versionlog in ?"bipartite-package". For citation see: citation("bipartite").
##  Have a nice time plotting and analysing two-mode networks.

## 
## Attaching package: 'bipartite'

## The following object is masked from 'package:vegan':
## 
##     nullmodel

data(mosquin1967)
library(vegan)
library(phytools)

## Loading required package: ape

## 
## Attaching package: 'ape'

## The following objects are masked from 'package:sna':
## 
##     consensus, degree

## Loading required package: maps

## 
## Attaching package: 'phytools'

## The following object is masked from 'package:vegan':
## 
##     scores

library(taxize)
#remotes::install_github("ropensci/brranching")
library(brranching)
library(rotl)

## 
## Attaching package: 'rotl'

## The following objects are masked from 'package:taxize':
## 
##     synonyms, tax_name, tax_rank

library(dendextend)

## Registered S3 method overwritten by 'dendextend':
##   method     from 
##   rev.hclust vegan

## 
## ---------------------
## Welcome to dendextend version 1.17.1
## Type citation('dendextend') for how to cite the package.
## 
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
## 
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## You may ask questions at stackoverflow, use the r and dendextend tags: 
##   https://stackoverflow.com/questions/tagged/dendextend
## 
##  To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
## ---------------------

## 
## Attaching package: 'dendextend'

## The following object is masked from 'package:phytools':
## 
##     untangle

## The following objects are masked from 'package:ape':
## 
##     ladderize, rotate

## The following object is masked from 'package:permute':
## 
##     shuffle

## The following object is masked from 'package:data.table':
## 
##     set

## The following object is masked from 'package:stats':
## 
##     cutree

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.0     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between()     masks data.table::between()
## ✖ dplyr::filter()      masks stats::filter()
## ✖ dplyr::first()       masks data.table::first()
## ✖ lubridate::hour()    masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag()         masks stats::lag()
## ✖ dplyr::last()        masks data.table::last()
## ✖ purrr::map()         masks maps::map()
## ✖ lubridate::mday()    masks data.table::mday()
## ✖ lubridate::minute()  masks data.table::minute()
## ✖ lubridate::month()   masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second()  masks data.table::second()
## ✖ purrr::transpose()   masks data.table::transpose()
## ✖ lubridate::wday()    masks data.table::wday()
## ✖ lubridate::week()    masks data.table::week()
## ✖ dplyr::where()       masks ape::where()
## ✖ lubridate::yday()    masks data.table::yday()
## ✖ lubridate::year()    masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(heatmaply)

## Loading required package: plotly
## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
## 
## Loading required package: viridis
## Loading required package: viridisLite
## 
## Attaching package: 'viridis'
## 
## The following object is masked from 'package:maps':
## 
##     unemp
## 
## Registered S3 method overwritten by 'seriation':
##   method         from 
##   reorder.hclust vegan
## 
## ======================
## Welcome to heatmaply version 1.5.0
## 
## Type citation('heatmaply') for how to cite the package.
## Type ?heatmaply for the main documentation.
## 
## The github page is: https://github.com/talgalili/heatmaply/
## Please submit your suggestions and bug-reports at: https://github.com/talgalili/heatmaply/issues
## You may ask questions at stackoverflow, use the r and heatmaply tags: 
##   https://stackoverflow.com/questions/tagged/heatmaply
## ======================

library(googlesheets4)

#load network

#load network
setwd("~/Desktop/Scent-Net")
gs4_auth(email=TRUE)
net.long <- read_sheet("1plhxbVh5IQLNVMazqO4OvfM3X6vbll_Y6BTMJ94bkI4", na=c("","na", "NA")) %>%    
      mutate(pollinator_group= if_else(pollinator=="papilio_gothica", "butterfly", pollinator_group)) %>% #fix papilio_gothica from bee -> butterfly
      mutate(plant= recode(plant,senecio_tall_int = "senecio_integerrimus", agoceris_glauca = "agoseris_glauca", agoceris_aurantiaca="agoseris_aurantiaca", hydrophyllum_fendeleri="hydrophyllum_fendleri"))



net <- net.long %>% select(plant, pollinator, interactions) %>% pivot_wider(names_from = plant, values_from = interactions, values_fn= sum, values_fill= 0) %>%  #pivots data set longer so plants as columns and pollinators as rows, replaces NAs with zeros
    column_to_rownames("pollinator") #pollinators are in first column, makes them into rownames; #deletes the first pollinators column


net <- as.matrix(t(net)) # turns it from dataframe to matrix bc packages like matrix; t transposes (turns dataset sideways - rows -> columns; columns -> rows) 
polls <- colnames(net) #now pollinators are column names
plants<-rownames(net)
snet <- sortweb(net) #sort by row & column totals; sorted rows by rowsums and columns by column sums. Biggest numbers in top left of matrix

#load scent

#Go to comm_vols for code to make scent.net
#Loads in vol, filters with bouquet, takes mean of each voc for each species 
#From comm_vols.Rmd
scent.net <- read_csv("data/Volatiles_by_species_mean.csv") %>% column_to_rownames("species") %>% #species are now rownames
    as.matrix()   #network analysis likes matrix 
    
scent.net <- scent.net[,colnames(scent.net)!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-"] #filtering this compound bc it has 0 interactions or emissions... 355 vs 354. breaks code later

#loads in chemical class and gets rid of ?
compound.class <-read_sheet("1BlDIOs4STe5ZTTPYE6qez1Nsey_ovEpe5m-SN8Q55gg", sheet= "chemical_class") %>% mutate(major_class= factor(major_class,levels=c("aliphatic", "monoterpene", "sesquiterpene", "nitrogenous", "sulfur", "benzenoid")) %>% fct_na_value_to_level("other")) %>% 
      filter(compound_long!="Cyclohexane, 1,3,5-trimethyl-2-octadecyl-") #filtering this compound bc it has 0 interactions or emissions... 355 vs 354

#codes chemical class by colors
voc.class.colors <- setNames(brewer.pal(7,"Set3"), levels(compound.class$major_class))

##match paul’s plants with the plants we have vocs for

scentplants<-rownames(scent.net) #%>% tolower() %>%  str_replace(" ", "_")

paulsplants <- plants

#plants not in both datasets 
setdiff(scentplants, paulsplants) #gives us bonus plants

   [1] "aquilegia_caerulea" "frasera_speciosa"   "geum_aleppicum"

setdiff(paulsplants, scentplants) 

   [1] "hydrophyllum_fendleri" "phacelia_heterophylla"

bothplants <- intersect(scentplants, paulsplants) #plants found in amanda and paul data

Pollinator groups and family

polls.count <-count(net.long, pollinator_group,pollinator_family, pollinator) #counting number of rows for each pollinator
polls.class <- tibble(pollinator=polls) %>%  left_join(polls.count) %>% pull(pollinator_group)
polls.class.dataframe <- tibble(pollinator=polls) %>%  left_join(polls.count)
polls.family <- tibble(pollinator=polls) %>%  left_join(polls.count) %>% pull(pollinator_family)
polls.order.all <- polls.class

#net.grp.all has all of pauls 45 plant species 
#net.grp only has plant species found in Amanda's and paul's datasets
net.grp.all <- aggregate(.~ polls.class, as.data.frame(t(net)), sum) #lump network by pollinator class; plants columns, pollinators rows
rownames(net.grp.all) <- net.grp.all[,1] #make plants rows, pollinators column names
net.grp.all[,1] <- NULL
net.grp.all <- as.matrix(t(net.grp.all)) #transposes matrix
net.grp<- net.grp.all[bothplants,] 

net.cut <- net[bothplants,]

Plant-pollinator network

##Pollinator groups

#Tells you plants major pollinator
plants.majorpoll <- setNames(factor(colnames(net.grp)[apply(net.grp, 1, which.max)]), rownames(net.grp))

#sets colors for ALL 5 pollinator groups
pcols.all <- setNames(brewer.pal(5,"Set2"), colnames(net.grp))

#Four colors for each major pollinator group
pcols <- pcols.all[ levels(plants.majorpoll)]

Proportion of visits by pollinator group

#gives % bee pollination for each plant species 
plants.bees <- net.grp[,"bee"]/rowSums(net.grp)
plot(sort(plants.bees))

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for

# Cap of plants.bees
plant.bee.cap <- capscale(sqrt(scent.net.cut)~plants.bees)
anova(plant.bee.cap)

   Permutation test for capscale under reduced model
   Permutation: free
   Number of permutations: 999
   
   Model: capscale(formula = sqrt(scent.net.cut) ~ plants.bees)
            Df Variance      F Pr(>F)
   Model     1   505636 1.0525   0.35
   Residual 41 19696248

#p=0.348
plot(plant.bee.cap, type="text")

#gives % fly pollination for each plant species 
plants.fly <- net.grp[,"fly"]/rowSums(net.grp)
plot(sort(plants.fly))

sort(plants.fly)

          agoseris_aurantiaca      amelanchier_alnifolia 
                  0.000000000                0.000000000 
         castilleja_sulphurea            gentiana_parryi 
                  0.000000000                0.000000000 
          ipomopsis_aggregata       lathyrus_lanszwertii 
                  0.000000000                0.000000000 
               mahonia_repens          mertensia_ciliata 
                  0.000000000                0.000000000 
         ranunculus_inamoenus               rosa_woodsii 
                  0.000000000                0.000000000 
              vicia_americana            viola_praemorsa 
                  0.000000000                0.000000000 
      delphinium_nuttallianum     hydrophyllum_capitatum 
                  0.001774623                0.010204082 
       campanula_rotundifolia       mertensia_fusiformis 
                  0.011428571                0.011627907 
     erythronium_grandiflorum      heliomeris_multiflora 
                  0.019230769                0.033676658 
            sedum_lanceolatum        heterotheca_villosa 
                  0.035398230                0.044919786 
           lomatium_dissectum       taraxacum_officinale 
                  0.076923077                0.102841678 
                senecio_serra       senecio_integerrimus 
                  0.125000000                0.174044876 
            arenaria_congesta           boechera_stricta 
                  0.178947368                0.181818182 
           erigeron_speciosus      solidago_multiradiata 
                  0.202226345                0.205357143 
         claytonia_lanceolata helianthella_quinquenervis 
                  0.244186047                0.267904509 
          fragaria_virginiana              linum_lewisii 
                  0.333333333                0.333333333 
          potentilla_hippiana       erigeron_flagellaris 
                  0.397425583                0.432900433 
          potentilla_gracilis       eriogonum_umbellatum 
                  0.433962264                0.442477876 
    pseudocymopterus_montanus                draba_aurea 
                  0.490196078                0.564285714 
         eriogonum_subalpinum  androsace_septentrionalis 
                  0.575000000                0.630000000 
               galium_boreale       achillea_millefolium 
                  0.875000000                0.942622951 
              agoseris_glauca 
                  1.000000000

#gives % butterfly pollination for each plant species 
plants.butterfly <- net.grp[,"butterfly"]/rowSums(net.grp)
plot(sort(plants.butterfly))

sort(plants.butterfly)

              agoseris_glauca      amelanchier_alnifolia 
                  0.000000000                0.000000000 
    androsace_septentrionalis     campanula_rotundifolia 
                  0.000000000                0.000000000 
         castilleja_sulphurea                draba_aurea 
                  0.000000000                0.000000000 
     erythronium_grandiflorum        fragaria_virginiana 
                  0.000000000                0.000000000 
              gentiana_parryi     hydrophyllum_capitatum 
                  0.000000000                0.000000000 
                linum_lewisii         lomatium_dissectum 
                  0.000000000                0.000000000 
               mahonia_repens          mertensia_ciliata 
                  0.000000000                0.000000000 
         mertensia_fusiformis  pseudocymopterus_montanus 
                  0.000000000                0.000000000 
         ranunculus_inamoenus               rosa_woodsii 
                  0.000000000                0.000000000 
            sedum_lanceolatum              senecio_serra 
                  0.000000000                0.000000000 
              vicia_americana            viola_praemorsa 
                  0.000000000                0.000000000 
        heliomeris_multiflora        potentilla_hippiana 
                  0.003934730                0.006436042 
      delphinium_nuttallianum       taraxacum_officinale 
                  0.006654836                0.008119080 
         achillea_millefolium helianthella_quinquenervis 
                  0.012295082                0.013262599 
        solidago_multiradiata        potentilla_gracilis 
                  0.013392857                0.018867925 
          heterotheca_villosa       claytonia_lanceolata 
                  0.022887701                0.023255814 
         lathyrus_lanszwertii          arenaria_congesta 
                  0.025000000                0.057142857 
         senecio_integerrimus       erigeron_flagellaris 
                  0.065494239                0.089466089 
         eriogonum_subalpinum             galium_boreale 
                  0.125000000                0.125000000 
         eriogonum_umbellatum         erigeron_speciosus 
                  0.150442478                0.179962894 
          ipomopsis_aggregata           boechera_stricta 
                  0.217391304                0.363636364 
          agoseris_aurantiaca 
                  1.000000000

#gives % hummingbird pollination for each plant species 
plants.hummingbird <- net.grp[,"hummingbird"]/rowSums(net.grp)
plot(sort(plants.hummingbird))

sort(plants.hummingbird)

         achillea_millefolium        agoseris_aurantiaca 
                  0.000000000                0.000000000 
              agoseris_glauca      amelanchier_alnifolia 
                  0.000000000                0.000000000 
    androsace_septentrionalis          arenaria_congesta 
                  0.000000000                0.000000000 
       campanula_rotundifolia       castilleja_sulphurea 
                  0.000000000                0.000000000 
         claytonia_lanceolata                draba_aurea 
                  0.000000000                0.000000000 
         erigeron_flagellaris         erigeron_speciosus 
                  0.000000000                0.000000000 
         eriogonum_subalpinum       eriogonum_umbellatum 
                  0.000000000                0.000000000 
     erythronium_grandiflorum        fragaria_virginiana 
                  0.000000000                0.000000000 
               galium_boreale            gentiana_parryi 
                  0.000000000                0.000000000 
   helianthella_quinquenervis      heliomeris_multiflora 
                  0.000000000                0.000000000 
          heterotheca_villosa     hydrophyllum_capitatum 
                  0.000000000                0.000000000 
         lathyrus_lanszwertii              linum_lewisii 
                  0.000000000                0.000000000 
           lomatium_dissectum             mahonia_repens 
                  0.000000000                0.000000000 
            mertensia_ciliata       mertensia_fusiformis 
                  0.000000000                0.000000000 
          potentilla_gracilis        potentilla_hippiana 
                  0.000000000                0.000000000 
    pseudocymopterus_montanus       ranunculus_inamoenus 
                  0.000000000                0.000000000 
                 rosa_woodsii          sedum_lanceolatum 
                  0.000000000                0.000000000 
                senecio_serra      solidago_multiradiata 
                  0.000000000                0.000000000 
         taraxacum_officinale            vicia_americana 
                  0.000000000                0.000000000 
              viola_praemorsa       senecio_integerrimus 
                  0.000000000                0.002425713 
             boechera_stricta    delphinium_nuttallianum 
                  0.090909091                0.224933452 
          ipomopsis_aggregata 
                  0.608695652

#gives % wasp pollination for each plant species 
plants.wasp <- net.grp[,"wasp"]/rowSums(net.grp)
plot(sort(plants.wasp))

sort(plants.wasp)

         achillea_millefolium        agoseris_aurantiaca 
                 0.0000000000               0.0000000000 
              agoseris_glauca      amelanchier_alnifolia 
                 0.0000000000               0.0000000000 
    androsace_septentrionalis          arenaria_congesta 
                 0.0000000000               0.0000000000 
             boechera_stricta     campanula_rotundifolia 
                 0.0000000000               0.0000000000 
         castilleja_sulphurea       claytonia_lanceolata 
                 0.0000000000               0.0000000000 
      delphinium_nuttallianum                draba_aurea 
                 0.0000000000               0.0000000000 
         erigeron_flagellaris         erigeron_speciosus 
                 0.0000000000               0.0000000000 
         eriogonum_umbellatum   erythronium_grandiflorum 
                 0.0000000000               0.0000000000 
          fragaria_virginiana             galium_boreale 
                 0.0000000000               0.0000000000 
              gentiana_parryi helianthella_quinquenervis 
                 0.0000000000               0.0000000000 
        heliomeris_multiflora     hydrophyllum_capitatum 
                 0.0000000000               0.0000000000 
          ipomopsis_aggregata       lathyrus_lanszwertii 
                 0.0000000000               0.0000000000 
                linum_lewisii         lomatium_dissectum 
                 0.0000000000               0.0000000000 
               mahonia_repens          mertensia_ciliata 
                 0.0000000000               0.0000000000 
         mertensia_fusiformis        potentilla_hippiana 
                 0.0000000000               0.0000000000 
    pseudocymopterus_montanus       ranunculus_inamoenus 
                 0.0000000000               0.0000000000 
                 rosa_woodsii          sedum_lanceolatum 
                 0.0000000000               0.0000000000 
         senecio_integerrimus              senecio_serra 
                 0.0000000000               0.0000000000 
         taraxacum_officinale            vicia_americana 
                 0.0000000000               0.0000000000 
              viola_praemorsa        heterotheca_villosa 
                 0.0000000000               0.0005347594 
        solidago_multiradiata        potentilla_gracilis 
                 0.0178571429               0.0377358491 
         eriogonum_subalpinum 
                 0.0750000000

##Heatmap1

#pp_heatmap
heatmaply(net, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          row_side_palette=function(n) pcols, 
          col_side_colors=data.frame(PollinatorGroup=polls.class), 
          col_side_palette=function(n) pcols.all)

##Matrix plot

#interaction matrix plot
visweb(net, type="nested", circles=T, boxes=F, circle.max=1.2)

##Web plot

#fix names of species (plants rows, polls cols)
net.nice.labels <- net
rownames(net.nice.labels) <- str_to_sentence(str_replace_all(rownames(net), "_"," "))
colnames(net.nice.labels) <- str_to_sentence(str_replace_all(colnames(net), "_"," "))

#interaction web plot
par(mar=c(0,0,0,0))
plant.pollinator.network <- plotweb(sqrt(net.nice.labels), empty=F, text.rot=90, method="cca", col.high=pcols.all[polls.order.all], col.low=pcols[plants.majorpoll], col.interaction=pcols.all[polls.order.all], bor.col.interaction=pcols.all[polls.order.all], bor.col.high=pcols.all[polls.order.all], bor.col.low=pcols[plants.majorpoll])

library(ggplot2)
ggsave(plant.pollinator.network, file='plant.pollinator.networkplot.jpg', width=130, height=105, units= 'mm', dpi=600)

Network topography

#computes network statistics (like nestedness)
#networkstats <- networklevel(net)
#write_csv(enframe(networkstats), "data/p-p network stats.csv")

#computes network statistics for voc-pollinator network
#scentnetworkstats<- networklevel(scentpollnet)
#write_csv(enframe(scentnetworkstats), "data/scent network stats.csv")
scentnet.stats <- read_csv("data/scent network stats.csv")

#the network statistics null model: TODO- broken
#network.nullmodel<-nullmodel(scentpollnet) 
#nullmodel.stats<- map(network.nullmodel, ~networklevel(.x )) #network.nullmodel[1:2]; index=c("NODF", "weighted NODF", "connectance", "H2", "modularity Q")

#nullmodel.stats %>% map_dfr(~enframe(.x), .id = "run") %>% write_csv("data/scent network null model statistics.csv")
nullmodel.stats <- read_csv("data/scent network null model statistics.csv")

#idea - pivot wider so network properties are the columns, rows are run number, then average by values in column -> pivot longer or store in new variable
#nullmodel.stats %>% pivot_wider(names_from = "Name", values_from = "Value") %>% column_to_rownames("stats") %>% as.matrix()

null.NODF <- read_csv("data/scent network null model statistics.csv") %>% filter(name=="NODF")
ggplot(null.NODF, aes(x=value))+ geom_histogram() + geom_vline(aes(xintercept=filter(scentnet.stats, name=="NODF") %>% pull(value)))

Scent network

##Heatmap2

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ] #cuts scent to match what paul has pollinator data for
plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)] #cuts Pauls data to match Amanda's

plants.bees.cut<- plants.bees[rownames(scent.net.cut)]

heatmaply(sqrt(scent.net.cut), scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          row_side_colors = data.frame(major.polls=plants.majorpoll.cut),
          row_side_palette=function(n) pcols )

          #col_side_colors=data.frame(PollinatorGroup=scent.net), #put compound class here
        #  col_side_palette=function(n) pcols.all)

Distance Matrix (for bray curtis distance)

Mantel test

#scent-plant network 
distance.matrix <- vegdist(scent.net.cut)
View(as.matrix(distance.matrix)) #tells you how different the smells are. If 0 there is no difference, scents are the same; 1 = no compounds in common

#plant-pollinator network
#how similar are the plants in terms of their pollinator array/composition
distance.matrix.pp <- vegdist(net.cut)
View(as.matrix(distance.matrix.pp))

#Mantel Test
#can do correlation, if scent similarity predicts pollinator --> mantel test
mantel(distance.matrix,distance.matrix.pp)

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix, ydis = distance.matrix.pp) 
   
   Mantel statistic r: 0.02341 
         Significance: 0.321 
   
   Upper quantiles of permutations (null model):
      90%    95%  97.5%    99% 
   0.0690 0.0835 0.0971 0.1197 
   Permutation: free
   Number of permutations: 999

#TODO: group by family then rerun vegdist

distance.matrix.pollgroup <- vegdist(net.grp)
mantel(distance.matrix, distance.matrix.pollgroup) 

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup) 
   
   Mantel statistic r: -0.06377 
         Significance: 0.892 
   
   Upper quantiles of permutations (null model):
      90%    95%  97.5%    99% 
   0.0670 0.0922 0.1132 0.1510 
   Permutation: free
   Number of permutations: 999

#can't predict which pollinator group will visit plant based on VOCS


#relative amounts of visits from pollinator groups 
distance.matrix.pollgroup.relative <- vegdist(decostand(net.grp, method="tot"))
mantel(distance.matrix, distance.matrix.pollgroup.relative) 

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix, ydis = distance.matrix.pollgroup.relative) 
   
   Mantel statistic r: -0.1527 
         Significance: 0.983 
   
   Upper quantiles of permutations (null model):
     90%   95% 97.5%   99% 
   0.106 0.142 0.167 0.206 
   Permutation: free
   Number of permutations: 999

#similarity of the relative amounts of scent of each plant
distance.matrix.relative <- vegdist(decostand(scent.net.cut, method="tot"))

#is similarity of the relative amounts of scent of each plant correlated with the similarity in relative visits from each pollinator group
mantel(distance.matrix.relative, distance.matrix.pollgroup.relative)

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pollgroup.relative) 
   
   Mantel statistic r: -0.1544 
         Significance: 0.982 
   
   Upper quantiles of permutations (null model):
     90%   95% 97.5%   99% 
   0.109 0.149 0.177 0.210 
   Permutation: free
   Number of permutations: 999

#No correlation between scent similarity and similarity of the proportion of visits by pollinators or pollinator group (p> 0.1)


#not grouped by pollinator class but still relative
distance.matrix.pp.relative <- vegdist(decostand(net.cut, method="tot"))
mantel(distance.matrix.relative, distance.matrix.pp.relative)

   
   Mantel statistic based on Pearson's product-moment correlation 
   
   Call:
   mantel(xdis = distance.matrix.relative, ydis = distance.matrix.pp.relative) 
   
   Mantel statistic r: 0.05562 
         Significance: 0.181 
   
   Upper quantiles of permutations (null model):
      90%    95%  97.5%    99% 
   0.0767 0.1018 0.1208 0.1549 
   Permutation: free
   Number of permutations: 999

Web plot VOCs-Plants

#interaction web plot
#plant-voc bipartite
plotweb(sqrt(scent.net), text.rot=90, method="cca", col.high=voc.class.colors[compound.class$major_class], bor.col.high=voc.class.colors[compound.class$major_class], col.low=pcols[plants.majorpoll], bor.col.low=pcols[plants.majorpoll], col.interaction=pal[5], bor.col.interaction=pal[5])

Plant-scent-pollinator network

NMDS

scent.net.cut <- scent.net[rownames(scent.net)%in% paulsplants, ]
scent.mds <- metaMDS(sqrt(scent.net.cut), autotransform = F, trace=F)
scent.op <- ordiplot(scent.mds, type="n")
points(scent.op, what="species", pch=3, col="grey50")
points(scent.op, what="sites", cex=1.5, pch=19, col=pcols[plants.majorpoll])
text(scent.op, what="sites", cex=0.8, col=pcols[plants.majorpoll])

#are scents of plants with different major pollinators different? no... 
#do phylogeny, plants in same genus will cluster together need to account for phylogeny

cap

scent.cap <- capscale(sqrt(scent.net.cut)~plants.majorpoll.cut)
anova(scent.cap)

   Permutation test for capscale under reduced model
   Permutation: free
   Number of permutations: 999
   
   Model: capscale(formula = sqrt(scent.net.cut) ~ plants.majorpoll.cut)
            Df Variance      F Pr(>F)
   Model     3   808879 0.5422  0.812
   Residual 39 19393005

plot(scent.cap, type="text")

Indicator compounds

library(indicspecies)

plants.majorpoll.cut <- plants.majorpoll[rownames(scent.net.cut)]
  
mp <- multipatt(as.data.frame(scent.net.cut), plants.majorpoll.cut, control=how(nperm=999))
summary(mp)

   
    Multilevel pattern analysis
    ---------------------------
   
    Association function: IndVal.g
    Significance level (alpha): 0.05
   
    Total number of species: 354
    Selected number of species: 9 
    Number of species associated to 1 group: 4 
    Number of species associated to 2 groups: 3 
    Number of species associated to 3 groups: 2 
   
    List of species associated to each combination: 
   
    Group butterfly  #sps.  1 
                                                                   stat p.value  
   Carbamic acid, N-phenyl-, 1,5-dimethyl-1-vinyl-4-hexenyl ester 0.942   0.035 *
   
    Group hummingbird  #sps.  3 
                                                   stat p.value  
   Propanoic acid, 2-methyl-, 3-methylbutyl ester 1.000   0.043 *
   Petasitene                                     0.996   0.043 *
   2H-Pyran-2-one, tetrahydro-3,6-dimethyl-       0.962   0.037 *
   
    Group butterfly+hummingbird  #sps.  2 
                                                              stat p.value    
   Bicyclo[7.2.0]undec-4-ene, 4,11,11-trimethyl-8-methylene- 0.967   0.001 ***
   Methoxyacetic acid, tetradecyl ester                      0.941   0.019 *  
   
    Group fly+hummingbird  #sps.  1 
                         stat p.value  
   2-Butenal, 3-methyl- 0.921   0.045 *
   
    Group butterfly+fly+hummingbird  #sps.  2 
                          stat p.value   
   (E)-4-Oxohex-2-enal   0.941   0.008 **
   1-Propanol, 2-methyl- 0.890   0.030 * 
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Shannon diversity

###Do generalist plants have higher diversity of scent compounds?
plants.numlinks  <- rowSums(decostand(net.cut, "pa"))
#polls.classtree <- compute.brlen(polls.classtree$phylo,1)
#library(picante)
#plants.pdlinks   <- pd(net.lump, polls.classtree)$PD
plants.numcompounds <- rowSums(decostand(scent.net.cut, "pa"))
plants.sumcompounds <- rowSums(scent.net.cut)
plants.shannoncompounds <- diversity(scent.net.cut, index="shannon")
plants.shannonlinks <- diversity(net.grp, index="shannon")

scent.diversity <- tibble(numlinks=plants.numlinks, numcompounds=plants.numcompounds, sumcompounds=plants.sumcompounds, shannoncompounds=plants.shannoncompounds, shannonlinks=plants.shannonlinks)

ggplot(scent.diversity, aes(y=numcompounds, x=numlinks))+
    geom_point()+ geom_smooth()

anova(lm(numlinks~numcompounds, data=scent.diversity))

   Analysis of Variance Table
   
   Response: numlinks
                Df Sum Sq Mean Sq F value  Pr(>F)  
   numcompounds  1  792.3  792.30  3.8827 0.05556 .
   Residuals    41 8366.4  204.06                  
   ---
   Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

#tells us if plants with more diverse scent get more pollinator species
anova(lm(plants.numlinks~shannoncompounds, data=scent.diversity))

   Analysis of Variance Table
   
   Response: plants.numlinks
                    Df Sum Sq Mean Sq F value Pr(>F)
   shannoncompounds  1   35.4  35.438  0.1593 0.6919
   Residuals        41 9123.3 222.520

#plot(plants.numcompounds~plants.pdlinks)

#plot(plants.shannoncompounds~plants.pdlinks)

shannon.compoundxlinks.plot <- ggplot(scent.diversity, aes(y=shannonlinks, x=shannoncompounds))+
    geom_point()+ geom_smooth()
ggsave(shannon.compoundxlinks.plot, file='shannon.compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)

compoundxlinks.plot <- ggplot(scent.diversity, aes(y=plants.numlinks, x=shannoncompounds))+
    geom_point()+ geom_smooth()
ggsave(compoundxlinks.plot, file='compoundxlinks.plot.jpg', width=130, height=105, units= 'mm', dpi=300)

shannon.compound.plot <- ggplot(scent.diversity, aes(y=shannoncompounds, x=plants.majorpoll.cut))+
    geom_point()+ geom_smooth()
ggsave(shannon.compound.plot, file='shannon.compound.plot.jpg', width=110, height=105, units= 'mm', dpi=300)

numcompounds.plot <-ggplot(scent.diversity, aes(y=numcompounds, x=plants.majorpoll.cut))+
    geom_point()+ geom_smooth()
ggsave(numcompounds.plot, file='numcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)

sumcompounds.plot<- ggplot(scent.diversity, aes(y=sumcompounds, x=plants.majorpoll.cut))+
    geom_point()+ geom_smooth()
ggsave(sumcompounds.plot, file='sumcompounds.plot.jpg', width=110, height=105, units= 'mm', dpi=300)

VOC-pollinator network

#interaction web plot
#voc-pollinator bipartite

#make dataframe with vocs and pollinators
##change scent.net to long format
plantcompoundemission <- scent.net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
  pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)


#write_csv(plantcompoundemission, "plantcompoundemissions.csv")
 
#total number of visits to plant by each pollinator species
plantpollinatorvisits <- net.cut %>% as.data.frame() %>% rownames_to_column("plant") %>% pivot_longer(-plant, names_to = "pollinator", values_to = "visits")


#total number of visits to each compound. Pollinators rows, compounds cols, values = number of visits
scentpollnet <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>%  group_by(pollinator, compound) %>% summarise(visits=sum(visits)) %>%  pivot_wider(names_from = "compound", values_from = "visits") %>% column_to_rownames("pollinator") %>% as.matrix()


#relative amount of each compound
plantcompoundrelative <- scent.net.cut %>% as.data.frame() %>% decostand(method="tot") %>% rownames_to_column("plant") %>% #average peak area
  pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)

#multiply number of visits by relative peak area
# weights visit to compounds. Proportions visits based on how many VOCs there are.
scentpollnet.area <- left_join(plantcompoundrelative, plantpollinatorvisits,relationship = "many-to-many") %>%  group_by(pollinator, compound) %>% mutate(visits.area=visits*emissions) %>%  summarise(visits.area=sum(visits.area)) %>%  pivot_wider(names_from = "compound", values_from = "visits.area") %>% column_to_rownames("pollinator") %>% as.matrix()

scentpollnet.class <- tibble(pollinator=rownames(scentpollnet)) %>% left_join(polls.class.dataframe) %>% pull(pollinator_group)


#plot bipartite network
#plotweb(sqrt(scentpollnet), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], #col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], #bor.col.interaction=pal[5])

#plotweb for scentpollnet.area

plotweb(sqrt(scentpollnet.area), text.rot=90, method="cca", col.high=pal[2], bor.col.high=pal[2], col.low=pcols[scentpollnet.class], bor.col.low=pcols[scentpollnet.class], col.interaction=pal[5], bor.col.interaction=pal[5])

Heatmap scent-net

#value = number of 
heatmaply(scentpollnet, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          #row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          #row_side_palette=function(n) pcols, 
          row_side_colors=data.frame(PollinatorGroup=scentpollnet.class), 
          row_side_palette=function(n) pcols.all)

voc.pollgrp.cap <- capscale(sqrt(scentpollnet.area)~scentpollnet.class)
anova(voc.pollgrp.cap)

   Permutation test for capscale under reduced model
   Permutation: free
   Number of permutations: 999
   
   Model: capscale(formula = sqrt(scentpollnet.area) ~ scentpollnet.class)
            Df Variance    F Pr(>F)
   Model     4   15.904 1.67  0.111
   Residual 88  209.505

plot(voc.pollgrp.cap, type="text")

chemical class x pollinator group bipartite

#make vector with compound class and all scent compounds
vol.class <-deframe(select(compound.class, compound_long, major_class) )

#turns scentpollnet matrix to data frame then pivots longer to have 2 cols called chemical and visits

chemical.taxa <- left_join(plantcompoundemission, plantpollinatorvisits,relationship = "many-to-many") %>% left_join(polls.class.dataframe) %>% left_join(compound.class, by=c("compound"= "compound_long")) %>% group_by(pollinator_group, major_class) %>% summarise(visits=sum(visits)) %>% pivot_wider(names_from = major_class, values_from = visits) %>% column_to_rownames("pollinator_group") %>% as.matrix()



#chemical class & pollinator group network
# how do I incorporate chemical class instead of compounds? 
#make data frame with chemical class and pollinator family - use scent.net & chemical.class
#and use scentpollnet but scentpollnet compounds in different order than scent.net
plotweb(sqrt(chemical.taxa), text.rot=90, method="normal", col.high=pal[3], bor.col.high=pal[3], col.low=pcols.all, bor.col.low=pcols.all, col.interaction=pal[5], bor.col.interaction=pal[5])

library(bipartite)
library(vcd)

#side quest
mosaic(chemical.taxa)

mosaic(net)

chisq.test(chemical.taxa) #p=1 so this means that no specialist compound classes. and vise versa. nobody is specialized 

   
    Pearson's Chi-squared test
   
   data:  chemical.taxa
   X-squared = 8488.8, df = 24, p-value < 2.2e-16

relative peak area x visits

#value = number of 
heatmaply(scentpollnet.area, scale="none", showticklabels = c(F,T), margins= c(0,NA,0,0),
          scale_fill_gradient_fun = ggplot2::scale_fill_gradient(low = "white", high = "#000050", trans="log1p"),
          distfun=vegdist, 
          hclustfun=function(x) hclust(x, method="mcquitty"), 
          #row_side_colors = data.frame(P=plants.majorpoll[rowSums(net)>0]),
          #row_side_palette=function(n) pcols, 
          row_side_colors=data.frame(PollinatorGroup=scentpollnet.class), 
          row_side_palette=function(n) pcols.all)

Rarifaction

#species accumulation curves
#figure out what percent Pauls data catpured the whole pollinator community
#pull out 1 random 1 hr sampling period and say how many pollinators observered, then do this again and again randomly - variation in samples = standard error y axis= number of species, x= number of pollinator observations 

#class= bee, fly, butterfly, hummingbird
#xaxis is number of plant species observed (if you did the whole summer but only looked at 1 species and then next summer looked at 2 species, etc.)
#

#simplify net.long to only have numbers in it so that accumulation curve will work 
#rows are observation "units" and columns are pollinators
net.wide <- net.long %>% select(pollinator, interactions) %>% mutate(index=row_number()) %>% 
    pivot_wider(names_from=pollinator, values_from = interactions) %>% mutate(across(everything(),~replace_na(.x, 0))) %>% select(-index) %>% as.matrix()

for(class in unique(polls.class)) { #could replace class with family
  class_subset <- net[,polls.class==class] #making accumulation curve for each pollinator class
  class_sa <- vegan::specaccum(class_subset, permutations = 999) 
  plot.new()
  plot(class_sa, col = as.integer(factor(class))+1, add = class != "bee", lwd=2, 
       xlab="Pollinator Observations", ylab="Pollinator Species", main="Pollinator Species accumulation curve", sub=class)
}

legend("bottomright", unique(polls.class), 
       col = 2:5, pch = 19, title = "Class")

sa<- vegan::specaccum(net.wide,permutations=100, method="exact")
plot(sa) 

#Pauls total pollinator observations 8hrs/wk * 42wks = 336 hours total

visitationrate <- nrow(net.long)/336 #15.33 unique p-p observations per hr (observation can include multiple visits)
#observation=within a 15min period how many p-p interactions did we see, total number of links observed
2000/visitationrate # 130 hrs per treatment ( 2 treatments = 260 hrs per year)

   [1] 130.4854

#then based on my observation hours, what percent of the pollinator community could be captured

Normalized Degree

Normalized degree for p-p

library(coin)
library(rstatix)
ND_poll<-ND(net, normalised = T)$higher 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll.class <- tibble(ND=ND_poll, polls.class=polls.class) 

ND.poll.class$ND <- as.numeric(ND.poll.class$ND)
ND.poll.class$polls.class <- as.factor(ND.poll.class$polls.class)

kruskal_test(ND~polls.class, data=ND.poll.class)

   
    Asymptotic Kruskal-Wallis Test
   
   data:  ND by
     polls.class (bee, butterfly, fly, hummingbird, wasp)
   chi-squared = 4.0146, df = 4, p-value = 0.404

ggplot(ND.poll.class, aes(x= polls.class, y=ND)) +geom_boxplot()

Normalized degree for VOC-pollinator

library(rstatix)

#ND for pollinators in scent network
ND_poll_scent<-ND(scentpollnet, normalised = T)$lower 

#ND of VOC chemical class
ND_scent_poll<-ND(scentpollnet, normalised = T)$higher


#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
ND.poll <-  tibble(ND=ND_poll_scent, polls.class=polls.class)

#kruskal_test doesn't like characters
ND.poll$ND <- as.numeric(ND.poll$ND)
ND.poll$polls.class <- as.factor(ND.poll$polls.class)
  
  
kruskal_test(ND~polls.class, data=ND.poll)

   
    Asymptotic Kruskal-Wallis Test
   
   data:  ND by
     polls.class (bee, butterfly, fly, hummingbird, wasp)
   chi-squared = 2.6914, df = 4, p-value = 0.6107

ggplot(ND.poll, aes(x= polls.class, y=ND)) +geom_boxplot()

#creates dataframe with normalized degree for each VOC and chemical class (same as data.frame)
ND.scent <- tibble(ND=ND_scent_poll, vols.class=compound.class$major_class) 
  
kruskal_test(ND~vols.class, data=ND.scent)

   
    Asymptotic Kruskal-Wallis Test
   
   data:  ND by
     vols.class (aliphatic, monoterpene, sesquiterpene, nitrogenous, sulfur, benzenoid, other)
   chi-squared = 7.8313, df = 6, p-value = 0.2507

ggplot(ND.scent, aes(x= vols.class, y=ND)) +geom_boxplot()

Normalized visits - home cooked

#created a function similar to normalized degree(ND) but ND looks at presence/absence in the network where normalized visits (NV) looks at the normalized weighted visits/interactions. This accounts for variation in the amount each plant or compound is visited and give those with higher interaction frequency greater weight or contribution to the network  

normalized.visits <- function (web, normalised = TRUE) 
{

  websum <- sum(web) #takes sum of all interactions in the web
  dlower <- rowSums(web)
  dhigher <- colSums(web)
  Nlow <- Nhigh <- 2
  if (normalised) {
    Nlow <- websum 
    Nhigh <- websum
  }
  low <- dlower/Nlow #weighting the interactions by dividing rowsums by the sum of all links in the network
  names(low) <- rownames(web) #dividing colsums by the sum of all links in the network
  high <- dhigher/Nhigh
  names(high) <- colnames(web)
  list(lower = low, higher = high)
}

Weighted normalized visits for p-p (pollinators only)- home cooked

#Which pollinators have the most visits. of all the visits recorded, what percentage of them were by that pollinator (for a given plant what percentage of visits does it get out of the whole meadow, for a pollinator what percentage of visits does it give out of all visits)

#normalized visits for pollinators in scent net
NV_poll<-normalized.visits(net, normalised = T)$higher 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class <- tibble(NV=NV_poll, polls.class=polls.class) 

#kruskal_test won't work as.character, need to change to numeric and factor
NV.poll.class$NV <- as.numeric(NV.poll.class$NV)
NV.poll.class$polls.class <- as.factor(NV.poll.class$polls.class)

kruskal_test(NV~polls.class, data=NV.poll.class)

   
    Asymptotic Kruskal-Wallis Test
   
   data:  NV by
     polls.class (bee, butterfly, fly, hummingbird, wasp)
   chi-squared = 11.58, df = 4, p-value = 0.02076

ggplot(NV.poll.class, aes(x= polls.class, y=NV)) +geom_boxplot()

#only a few species making most of the visits
#1 bee responsible for over 40% of visits

Normalized visits for scent network

#normalized visits for pollinators in scent net
#which pollinator class visits the most amount of compounds weighted by the number of visits
NV_poll_scentnet<-normalized.visits(scentpollnet, normalised = T)$lower 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.poll.class.scentnet <- tibble(NV=NV_poll_scentnet, polls.class=polls.class) 
  
#kruskal_test won't work as.character, need to change to numeric and factor
NV.poll.class.scentnet$NV <- as.numeric(NV.poll.class.scentnet$NV)
NV.poll.class.scentnet$polls.class <- as.factor(NV.poll.class.scentnet$polls.class)

kruskal_test(NV~polls.class, data= NV.poll.class.scentnet)

   
    Asymptotic Kruskal-Wallis Test
   
   data:  NV by
     polls.class (bee, butterfly, fly, hummingbird, wasp)
   chi-squared = 2.9103, df = 4, p-value = 0.573

ggplot(NV.poll.class.scentnet, aes(x= polls.class, y=NV)) +geom_boxplot()

#super generalist bees

#normalized visits for VOCs in scent net
#which compound class are most visited
NV_voc<-normalized.visits(scentpollnet, normalised = T)$higher 

#creates dataframe with normalized degree for each pollinator and pollinator class (same as data.frame)
NV.vols.class <- tibble(NV=NV_voc, vols.class=compound.class$major_class) 
  

k_t.vols.class <- kruskal_test(NV~vols.class, data=NV.vols.class)


library(coin)
library(PMCMRplus) #for posthoc
library(PMCMR)

#posthoc_vols.class <- posthoc.kruskal.nemenyi.test(NV.vols.class$vols.class, NV.vols.class$NV, dist = "tukey")

ggplot(NV.vols.class, aes(x= vols.class, y=NV)) +geom_boxplot()

#tells you what vocs get the most visits
#no difference in visitation rate to compound classes

#Modularity

Modularity for plant-pollinator network

library(vcd)

modules.net <- computeModules(sqrt(net))

plotModuleWeb(modules.net, labsize=.5)

module.net.list <- listModuleInformation(modules.net)

flatten(module.net.list[[2]])

   [[1]]
    [1] "androsace_septentrionalis" "claytonia_lanceolata"     
    [3] "draba_aurea"               "erythronium_grandiflorum" 
    [5] "hydrophyllum_capitatum"    "lomatium_dissectum"       
    [7] "mertensia_fusiformis"      "pseudocymopterus_montanus"
    [9] "ranunculus_inamoenus"      "taraxacum_officinale"     
   [11] "viola_praemorsa"           "amelanchier_alnifolia"    
   [13] "fragaria_virginiana"       "hydrophyllum_fendleri"    
   [15] "linum_lewisii"            
   
   [[2]]
    [1] "andrena_cyanophila"     "lasioglossum_spp"       "sphecodes_sp_01"       
    [4] "syrphidae_spp"          "osmia_iridis"           "scathophagidae_sp_01"  
    [7] "halictus_confusus"      "halictus_virgatellus"   "halictidae_green_spp_1"
   [10] "eupeodes_volucris"      "andrena_sp_02"          "muscidae_sp_02"        
   [13] "syrphidae_02"           "rhamphomyia_sp_01"      "cartosyrphus_tarda"    
   [16] "callophrys_affinis"     "osmia_lignaria"         "cartosyrphus_sp_02"    
   [19] "tephritidae_01"         "polygonia_zephyrus"    
   
   [[3]]
   [1] "delphinium_nuttallianum" "ipomopsis_aggregata"    
   [3] "lathyrus_lanszwertii"    "solidago_multiradiata"  
   [5] "gentiana_parryi"         "vicia_americana"        
   
   [[4]]
   [1] "bombus_appositus"        "selasphorus_platycercus"
   [3] "coelioxus_funeraria"     "anthidium_emarginatum"  
   [5] "bombus_nevadensis"       "papilio_gothica"        
   [7] "sphex_spp_1"            
   
   [[5]]
   [1] "agoseris_aurantiaca"        "campanula_rotundifolia"    
   [3] "erigeron_speciosus"         "helianthella_quinquenervis"
   [5] "heliomeris_multiflora"      "heterotheca_villosa"       
   [7] "senecio_serra"              "castilleja_sulphurea"      
   [9] "phacelia_heterophylla"     
   
   [[6]]
    [1] "nymphalidae_spp"         "bombus_bifarius"        
    [3] "megachile_relativa"      "tachinidae_sp_01"       
    [5] "bombus_sylvicola"        "lycaenidae_spp"         
    [7] "anthophora_terminalis"   "bombus_flavifrons"      
    [9] "andrena_costillensis"    "bombus_californicus"    
   [11] "bombus_occidentalis"     "megachile_pugnata"      
   [13] "arctophila_flagrans"     "osmia_grindeliae"       
   [15] "pieridae_spp"            "speyeria_mormonia"      
   [17] "chrysotoxum_ventricosum" "megachile_melanophea"   
   [19] "psithyrus_insularis"     "lycaena_rubidus_sirius" 
   [21] "colias_alexandra"        "bombus_frigidus"        
   [23] "mesebrina_latreillei"    "syrphidae_spp_2"        
   [25] "syrphidae_spp_4"         "bombus_balteatus"       
   [27] "syrphidae_spp_1"         "phormia_spp_1"          
   [29] "megachile_frigida"       "lycaena_cupreus"        
   [31] "bombus_mixtus"           "osmia_tristela"         
   [33] "melissodes_confusa"     
   
   [[7]]
    [1] "agoseris_glauca"      "arenaria_congesta"    "eriogonum_subalpinum"
    [4] "eriogonum_umbellatum" "mahonia_repens"       "mertensia_ciliata"   
    [7] "potentilla_hippiana"  "rosa_woodsii"         "sedum_lanceolatum"   
   [10] "achillea_millefolium" "potentilla_gracilis" 
   
   [[8]]
    [1] "muscidae_spp_01"          "colletes_kincaidii"      
    [3] "halictus_rubicundus"      "andrena_transnigra"      
    [5] "halictidae_spp_01"        "panurginus_ineptus"      
    [7] "thricops_septentrionalis" "osmia_bucephala"         
    [9] "chrysis_coerulans"        "villa_eumenes"           
   [11] "glaucopsyche_lygdamus"    "pieris_mcdunnoughi"      
   [13] "andrena_vicinoides"       "colletes_nigrifrons"     
   [15] "melanstoma_caerulescens" 
   
   [[9]]
   [1] "boechera_stricta"     "erigeron_flagellaris" "galium_boreale"      
   [4] "senecio_integerrimus"
   
   [[10]]
    [1] "ochlodes_sylvanoides"  "systoechus_sp"         "osmia_proxima"        
    [4] "pieris_rapa"           "hoplitus_fulgida"      "osmia_subaustralis"   
    [7] "osmia_montana"         "osmia_coloradensis"    "sphaerophoria_robusta"
   [10] "melanostoma_kelloggi"  "syrphidae_spp_3"       "eupeodes_lapponicus"  
   [13] "euphydras_anicia"      "eristalis_latifrons"   "osmia_sp_sgo"         
   [16] "cryptopogon_sp_01"     "hoplitus_robusta"      "coenonympha_ochracea"

pollinator.modules <- tibble(module=1:5) %>% mutate(pollinator=map(module, ~module.net.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()

chisqr.polls <- chisq.test(pollinator.modules)
residules.chisqr.polls <- residuals(chisqr.polls)

mosaic(t(pollinator.modules), shade = T, gp=shading_hcl)

#test to see p-p network modules
#mosaic(net, shade = T, gp=shading_hcl)

Modularity for scent-pollinator network

modules.scentnet <- computeModules(sqrt(t(scentpollnet)))

#plotModuleWeb(modules.scentnet)

module.scentnet.list <- listModuleInformation(modules.scentnet)


pollinator.modules <- tibble(module=1:4) %>% mutate(pollinator=map(module, ~module.scentnet.list[[2]][[.x]][[2]])) %>% unnest(pollinator) %>% left_join(polls.class.dataframe) %>% count(module, pollinator_group) %>% pivot_wider(names_from = pollinator_group, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()

chisqr.polls <- chisq.test(pollinator.modules)
residules.chisqr.polls <- residuals(chisqr.polls)

mosaic(t(pollinator.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))

#chemical class modularity 


chemical.modules <- tibble(module=1:4) %>% mutate(compound_long=map(module, ~module.scentnet.list[[2]][[.x]][[1]])) %>% unnest(compound_long) %>% left_join(compound.class) %>% count(module, major_class) %>% pivot_wider(names_from = major_class, values_from = n, values_fill = 0) %>% column_to_rownames("module") %>% as.matrix()

mosaic(t(chemical.modules), shade = T, gp=shading_hcl, labeling = vcd::labeling_border(rot_labels = c(0, 0)))

Stacked bar chart - who’s visiting who

Plot 1 -pollinator perspective

#Stacked bar chart with flowers visited by each pollinator #Looks at the top 20 plants and top 75 pollinators (interactions square rooted to adjust scale on the figure) #x= interactions, y=pollinators, color=plant

Plot 2 - plants perspective

#Stacked bar chart with what pollinators visited each plant species #45 plant species, 18 pollinators #x= interactions, y=plant, color=pollinator

library(ggplot2)
library(tidyverse)
library(RColorBrewer)
mypalette= c(brewer.pal(8, "Set2"), brewer.pal(12, "Set3"))

#looks at top 20 plants and top 75 pollinators (interactions sqrt to adjust scale on figure)
#pollinator perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 18), pollinator=fct_lump_n(pollinator, 75)), aes(x=sqrt(interactions), y=pollinator, fill=plant)) + geom_col() + scale_fill_manual(values=mypalette)

#plant perspective
ggplot(net.long %>% mutate(plant=fct_lump_n(plant, 45), pollinator=fct_lump_n(pollinator, 18)), aes(x=sqrt(interactions), y=plant, fill=pollinator)) + geom_col() + scale_fill_manual(values=mypalette)

Phenology

#Uses Paul’s phenology data to create stackplot

Plot 1 - Flower Phenology

#Panels are years, x= week number, y= flower count, color= plant species ## Plot 2 - flower phenology heatmap ## Plot 3 - Pollinator Phenology #Panels are years, x= week number, y= interactions, color= pollinator ## Plot 4 - flower phenology heatmap ## Plot 5 - Scent Phenology #weighted average of flower scents and mean is weighted by how many flowers there are that week #shows relative amount of that compound out of the total scent in the meadow that week #each color is a compound, value= total emissions ## Plot 6 - heatmap scent phenology

phenology <- read_sheet("15Jnx8OWMzpakC0YvJWmZMmaRSnA0X-7HBjymDsvUfjY")

#stack plots and heatmaps 

#Flowers#
ggplot(phenology, aes(x=week_num, y=flower_count, fill= plant )) + geom_col() + scale_fill_manual(values=sample(rainbow(46))) +facet_wrap(vars(year))

#heatmap
ggplot(phenology, aes(x=week_num, y=plant, fill= sqrt(flower_count ))) + geom_tile() + scale_fill_viridis_c()  + facet_wrap(vars(year))

#pollinators#

ggplot(net.long, aes(x=week_num, y=interactions, fill= pollinator )) + geom_col() + scale_fill_manual(values=sample(rainbow(93))) +facet_wrap(vars(year))

#heatmap
ggplot(net.long, aes(x=week_num, y=pollinator, fill= sqrt(interactions ))) + geom_tile() + scale_fill_viridis_c()  + facet_wrap(vars(year))

#scent#
#number of flowers per week and multiple by mean scent of that species

#how many flowers are open each week
flowers.per.week <- phenology %>% group_by(year, week_num, plant) %>%  
                    summarise(flower_count=sum(flower_count))

#ideally want for each week, 1 average weighted mean of scent compound 
  #but we only have a few scent samples so just taking the mean of scent compounds for each plant

#flowers per week percentage 
flower.percentage <- flowers.per.week %>% group_by(year, week_num) %>% 
        mutate(flower_count=flower_count/sum(flower_count)) %>% 
        pivot_wider(names_from = c(year,week_num), values_from = flower_count, values_fill = 0) %>% 
        filter(plant%in%rownames(scent.net)) %>% 
        arrange(match(plant, rownames(scent.net)))%>%  #matches scent.net plant order
        column_to_rownames("plant")  %>% as.matrix()
#change names in pauls data

scentxphenology <- t(flower.percentage) %*% scent.net[rownames(flower.percentage),]  
#values in matrix = weighted average of flower scents and mean is weighted by how many flowers there are that week

#replaces Nans with 0
scentxphenology[is.nan(scentxphenology)] <- 0

#ggplot likes long format
scent.phenology <- scentxphenology %>% as.data.frame() %>%  rownames_to_column("yearweek") %>% separate(yearweek, into = c("year", "week_num")) %>% mutate(year=as.integer(year), week_num= as.integer(week_num)) %>%  pivot_longer(all_of(colnames(scent.net)))

ggplot(scent.phenology, aes(x=week_num, y=value, fill= name )) + geom_col() + scale_fill_manual(values=sample(rainbow(355)), guide="none") +facet_wrap(vars(year))

#heatmap - scent x phenology
ggplot(scent.phenology, aes(x=week_num, y=name, fill= sqrt(value ))) + geom_tile() + scale_fill_viridis_c()  + facet_wrap(vars(year))

Mega table

flower and poll traits

VOCS

#number of compounds each plant species has
#total emissions for each species
#diversity of scents (shannon diversity)
library(tidyverse)

scent.net.long <- scent.net %>% as.data.frame() %>% rownames_to_column("plant") %>% #average peak area
  pivot_longer(-plant, names_to = "compound", values_to = "emissions") %>% filter(emissions!=0)


compound.class.percentage <- scent.net.long %>% left_join(compound.class, by=c(compound="compound_long")) %>%  group_by(plant, major_class) %>%  summarise(emissions=sum(emissions)) %>%  #gives total emissions for each compound class
        mutate(emissions=emissions/sum(emissions)) %>%  #relative amount of each compound class in species scent
        pivot_wider(names_from = major_class, values_from = emissions, names_prefix = "percent_", values_fill = 0) #makes column for each compound class and gives percentage of that compound class out of the total emissions for that species 


#rare vs common compounds (mean rarity)
##for each compound, how many times does it occur in the dataset? 

compound.rarity <- scent.net.long %>%  count(compound) %>%  #number of species the compound occurs in
          mutate(rarity=n/nrow(scent.net))  #the percentage of species the compound occurs in
 
weighted.scent.rarity <- scent.net.long %>% left_join(compound.rarity) %>% group_by(plant) %>%    
            mutate(emissions=emissions/sum(emissions)) %>%  #makes emissions -> relative emissions    
            mutate(weighted_scent_rarity = emissions*rarity) %>% #multiplies emissions x rarity to get weighted average
            summarise(weighted_scent_rarity= sum(weighted_scent_rarity)) #then add the different parts of the average
 #rarity score of 1 means this compound is found in all the plant species, 0= no other species has these compounds 
#mertensia has 0.70 rarity - the weighted average of all the compounds rarity; 70% of compounds (or total emissions) are shared, 30% of compounds or total emissions other species don't have



scent.summary <- tibble(plant=rownames(scent.net), total_emissions= rowSums(scent.net), number_compounds=rowSums(scent.net>0),scent_shannon_diversity= diversity(scent.net, index="shannon")) %>%
  left_join(weighted.scent.rarity) %>%     
  left_join(compound.class.percentage) 

floral traits summary

#pollinators that visit plants

#Percentage of visits by each pollinator family to each plant species
pollinator.family.percentage <- net.long %>% group_by(plant, pollinator_family) %>%  summarise(interactions=sum(interactions)) %>%  #gives sum of interactions for each pollinator family
        mutate(interactions=interactions/sum(interactions)) %>%  #relative amount of each pollinator family visiting each flower species
        pivot_wider(names_from = pollinator_family, values_from = interactions, names_prefix = "percent_", values_fill = 0)

#Percentage of visits by each pollinator group to each plant species
pollinator.group.percentage <- net.long %>% group_by(plant, pollinator_group) %>%  summarise(interactions=sum(interactions)) %>%  #gives sum of interactions for each pollinator group
        mutate(interactions=interactions/sum(interactions)) %>%  #relative amount of each pollinator group visiting each flower species
        pivot_wider(names_from = pollinator_group, values_from = interactions, names_prefix = "percent_", values_fill = 0)

#rarity - how generalized the pollinator is
pollinator.rarity <- net.long %>% count(plant, pollinator) %>%  #adds how many rows of interactions between plant & poll species
          count(pollinator, name="number_plants_visited") %>%  #number of plants the pollinator goes to
          mutate(rarity=number_plants_visited/nrow(net))  #the percentage of plant species the pollinator visits in the meadow

#gives for each plant species, which pollinators visited and their relative contributions to total number of visits
 net.long.relative.interactions <- net.long %>% group_by(plant, pollinator) %>% 
            summarise(interactions=sum(interactions)) %>% 
            left_join(pollinator.rarity) %>% 
            group_by(plant) %>%    
            mutate(interactions=interactions/sum(interactions))
 
 #percent of interactions that come from each pollinator species    
 weighted.pollinator.rarity <-   net.long.relative.interactions %>% mutate(weighted_poll_rarity = interactions*rarity) %>% #multiplies interactions x rarity to get weighted average
            summarise(weighted_poll_rarity = sum(weighted_poll_rarity)) 

flower.visitor.summary <- tibble(plant=rownames(net), total_visits= rowSums(net), number_pollspecies_visited=rowSums(net>0),pollinator_shannon_diversity= diversity(net, index="shannon")) %>%
  left_join(weighted.pollinator.rarity) %>% 
  left_join(pollinator.family.percentage) %>% left_join(pollinator.group.percentage)

megatable <- scent.summary %>% full_join(flower.visitor.summary) #full join bc pauls data and scent data plants don't match

Correlations

Plot 1 - heatmap of correlation table

Plot 2 - weighted scent and weighted pollinator visits correlation scatter plot

Plot 3 - percent_hesperiidae x nitrogenous vocs correlation scatter plot

library(pheatmap)
library(smatr)


correlation.table <- megatable %>% column_to_rownames("plant") %>% cor(use="na.or.complete")

pheatmap(correlation.table, scale="none", clustering_method = "mcquitty", color=colorRampPalette(c("blue","white", "red"))(200),
         breaks = seq(-1,1, by=0.01)) #which ones it labels on the scale

#scatter plot with just weighted scent and weighted pollinator visits
weighted.scent.poll <- ggplot(megatable, aes(x=weighted_poll_rarity, y=weighted_scent_rarity, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))     

print(weighted.scent.poll)

#the most common volatiles get the most visits from the most generalized pollinators

sma.weighted.scent.poll <- sma(weighted_poll_rarity~weighted_scent_rarity, data=megatable )

summary(sma.weighted.scent.poll)

   Call: sma(formula = weighted_poll_rarity ~ weighted_scent_rarity, data = megatable) 
   
   Fit using Standardized Major Axis 
   
   ------------------------------------------------------------
   Coefficients:
                elevation    slope
   estimate    -0.3602556 1.352232
   lower limit -0.5885207 1.001742
   upper limit -0.1319906 1.825352
   
   H0 : variables uncorrelated
   R-squared : 0.06769348 
   P-value : 0.091989

#percent_hesperiidae x nitrogenous vocs
ggplot(megatable, aes(x=percent_nitrogenous, y=percent_lycaenidae, color=number_pollspecies_visited)) + geom_point()+geom_smooth() + geom_text(aes(label=plant))

##pollinator trait summary